Windows 95 Super Collection/Windows 95 Super Collection.iso/win95/bench/thread/fpcalc.c

C/C++ Source or Header | 1995-08-13 | 4.4 KB | 192 lines

//==========================================================================
// FPCALC.C - C routines that run a floating-point intensive operation.
//             (Change this code as desired to implement your own floating-
//                 point benchmark.)
//
//                    Note: This floating-point calculation (a Fast Fourier
//                            transformation is borrowed from the
//                          public domain Stanford UNIX benchmarks.)
//
//                            DoFloatingPointCalc() is the entry point here,
//                   which is called by each thread procedure.
//
//==========================================================================

#include <stdlib.h>
#include "fpcalc.h"

UINT DoFloatingPointCalc(UINT uReps, UINT uRepsThisPass, UINT uCurrentReps)
{
    // Here, FP calculation is a fast Fourier transformation
    //       using matrices of complex numbers.
    UINT i = 0;

    while(1)
    {
        Oscar();

         i++;
        if(i >= uRepsThisPass || i + uCurrentReps >= uReps)
            break;
    }

    return i + uCurrentReps;
}

// Fast Fourier transform define's
#define fftsize           256
#define fftsize2          129

//----------------------------------------------------------
//  Fast Fourier Transformation (Oscar)
//----------------------------------------------------------
//----------------------------------------------------------
// Cos() computes cos of x (x in radians) by an expansion.
//----------------------------------------------------------
float Cos (float x)
{
    int i, factor;
    float    result,power;

   result = (float)1.0; factor = 1;  power = x;
   for ( i = 2; i <= 10; i++ )
   {
      factor = factor * i;  power = power*x;
      if ( (i & 1) == 0 )
      {
             if ( (i & 3) == 0 ) result = result + power/factor;
             else result = result - power/factor;
         }
   }
   return (result);
}

//----------------------------------------------------------
// Uniform11()
//----------------------------------------------------------
void Uniform11(long iy, float yfl)
{
    iy = (4855*iy + 1731) & 8191;
    yfl = (float)iy/(float)8192.0;
}

//----------------------------------------------------------
// Exptab()
//----------------------------------------------------------
void Exptab(int n, struct complex e[])
{
    float theta, divisor, h[26];
   int i, j, k, l, m;

   theta = (float)3.1415926536;
   divisor = (float)4.0;
   for ( i=1; i <= 25; i++ )
   {
       h[i] = 1/(2*Cos( theta/divisor ));
      divisor = divisor + divisor;
   }

   m = n / 2 ;
   l = m / 2 ;
   j = 1 ;
   e[1].rp = (float)1.0 ;
   e[1].ip = (float)0.0;
   e[l+1].rp = (float)0.0;
   e[l+1].ip = (float)1.0 ;
   e[m+1].rp = (float)-1.0 ;
   e[m+1].ip = (float)0.0 ;

   do
   {
       i = l / 2 ;
      k = i ;
         do
         {
             e[k+1].rp = h[j]*(e[k+i+1].rp+e[k-i+1].rp) ;
             e[k+1].ip = h[j]*(e[k+i+1].ip+e[k-i+1].ip) ;
             k = k+l ;
         }
         while ( k <= m );

         j = min( j+1, 25);
         l = i ;
   }
   while ( l > 1 );

}

//----------------------------------------------------------
// Fft()
//----------------------------------------------------------
void Fft( int n, struct complex z[], struct complex w[], struct complex e[], float sqrinv)
{
    int i, j, k, l, m, index;
   m = n / 2 ;
   l = 1 ;

   do
   {
         k = 0 ;
         j = l ;
         i = 1 ;

        do
        {
            do
            {
             w[i+k].rp = z[i].rp+z[m+i].rp ;
             w[i+k].ip = z[i].ip+z[m+i].ip ;
             w[i+j].rp = e[k+1].rp*(z[i].rp-z[i+m].rp)
                 -e[k+1].ip*(z[i].ip-z[i+m].ip) ;
             w[i+j].ip = e[k+1].rp*(z[i].ip-z[i+m].ip)
                 +e[k+1].ip*(z[i].rp-z[i+m].rp) ;
                 i = i+1 ;
             }
             while ( i <= j );

           k = j ;
             j = k+l ;
         }
         while ( j <= m );

         index = 1;
         do
         {
             z[index] = w[index];
             index = index+1;
         }
         while ( index <= n );
         l = l+l ;
   }
   while ( l <= m );

     for ( i = 1; i <= n; i++ )
     {
         z[i].rp = sqrinv*z[i].rp ;
         z[i].ip = -sqrinv*z[i].ip;
     }
}

void Oscar(void)
{
    struct complex z[fftsize+1], w[fftsize+1], e[fftsize2+1];
    float    zr, zi;

    int i;
    long seed = 5767 ;

    Exptab(fftsize,e) ;
    for ( i = 1; i <= fftsize; i++ )
    {
       Uniform11( seed, zr );
        Uniform11( seed, zi );
        z[i].rp = (float)20.0*zr - (float)10.0;
        z[i].ip = (float)20.0*zi - (float)10.0;
   }

    for ( i = 1; i <= 20; i++ )
    {
      Fft(fftsize,z,w,e,(float)0.0625) ;
      }
}